Effects of Thermal Noise on Pattern Onset in Continuum Simulations of Shaken 
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The author investigates the onset of patterns in vertically oscillated layers of dissipative particles 
using numerical solutions of continuum equations to Navier-Stokes order. Above a critical acceler- 
ational amplitude of the cell, standing waves form stripe patterns which oscillate subharmonically 
with respect to the cell. Continuum simulations neglecting interparticle friction yield pattern wave- 
lengths consistent with experiments using frictional particles. However, the critical acceleration for 
standing wave formation is approximately 10% lower in continuum simulations without added noise 
than in molecular dynamics simulations. This report incorporates fluctuating hydrodynamics theory 
into continuum simulations by adding noise terms with no fit parameters; this modification yields a 
critical acceleration in agreement with molecular dynamics simulations. 
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A successful theory of granular hydrodynamics would 
allow scientists and engineers to apply the powerful meth- 
ods of fluid dynamics to granular flow. Despite experi- 
mental [l], [1] and computational [1, 0] evidence demon- 
strating the potential utility of hydrodynamics models 
for grains, a general set of hydrodynamic governing equa- 
tions is not yet recognized for granular media [5|-[7|. 

One granular hydrodynamics approach derives contin- 
uum equations for number density n, velocity u, and 
granular temperature T (|T is the average kinetic en- 
ergy due to random particle motion) by modeling par- 
ticle interactions with binary, hard sphere collision op- 
erators in kinetic theory p4lO|- These equations repre- 
sent a different approach from other popular methods of 
modeling grains, such as molecular dynamics (MD) sim- 
ulations which simulate individual grain motion. This 
report is the first to directly incorporate fluctuating hy- 
drodynamics theory into continuum simulations of three- 
dimensional (3D) time-dependent granular flow. 

Vertically shaken layers provide an important testbed 
for granular phenomena |lll4l5l |. A flat layer of grains 
with depth H oscillated sinusoidally in the direction of 
gravity with frequency / and amplitude A leaves the 
plate at some time during the cycle if the maximum accel- 
eration of the plate amax — ^ (27r/)^ is greater than the 
acceleration of gravity g. Thus the layer leaves the plate 
if the dimensionless accelerational amplitude T = amax/g 
exceeds unity. When F exceeds a critical value Tc, the 
layer spontaneously forms standing waves which are sub- 
harmonic with respect to the plate. Various standing 
wave patterns are found experimentally, depending on F 
and the dimensionless frequency /* = f\/H/g ,14] . 

Previous experiments [l^ and MD simulations (T7| 
have shown that friction between grains plays a role in 
these patterns. Experimentally, adding graphite to re- 
duce friction decreased Fc and prevented the formation 
of stable square or hexagonal patterns found for certain 
ranges of /* and F in experiments without graphite jT6| . 
Similarly, MD simulations with friction between parti- 
cles have quantitatively reproduced stripe, square, and 



hexagonal subharmonic standing waves seen experimen- 
tally [1^ , but MD simulations without frictionyield only 
stable stripe patterns and display a lower Fc [17l | . In this 
report, I investigate the onset of stripe patterns in con- 
tinuum simulations of frictionless particles. 

Continuum equations for granular media have been 
proposed using a variety of approximations [l|, |^, 
I use a continuum simulation previously used to model 
shock waves [l9( and patterns in a granular shaker 
in order to directly compare to previous results The 
granular fluid is contained between two impenetrable hor- 
izontal plates at the top and bottom of the container. 
The lower plate oscillates sinusoidally between height 
z — and z = 2A, and the ceiling is located at a height 
Lz above the lower plate. Periodic boundary conditions 
are used in the horizontal directions x and y to eliminate 
sidewall effects. The dimensions of the box L^;, Ly, and 
Lz can be varied. This simulation numerically integrates 
equations of Navier-Stokes order proposed by Jenkins 
and Richman ^ for a dense gas of frictionless (smooth) , 
inelastic hard spheres with uniform diameter a. Energy 
loss due to collisions is characterized by a single param- 
eter, the normal coefficient of restitution e = 0.70. In- 
tegrating these hydrodynamic equations using a second 
order finite difference scheme on a uniform grid in 3D 
with first order adaptive time stepping [l^ yields num- 
ber density, momentum, and granular temperature. 

Above Fc, stripes are seen experimentally for a range 
of parameters, including nondimensional frequency /* = 
0.4174, and layer depth LL = 5.4fT [ll|. In this report, I 
compare to previous continuum and MD simulations [4\ , 
where F was varied while frequency /* — 0.4174 and the 
number of particles particles per unit area which 

experimentally corresponds to a layer depth LL = 5Aa as 
poured |18i] ) were fixed. This corresponds to a frequency 
of 56 Hz for particles with diameter a = 0.1mm. To com- 
pare current results to that previous investigation, I use 
the same frequency, layer depth, and cell size horizontally 
Lx = Ly ~ 42<T and vertically Lz — 80a [4|. 

In that report, continuum simulations produced flat 
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layers for accelerational amplitudes below r^"* = 
1.955±0.005, and stripe patterns above this critical value. 
MD simulations produced disordered peaks and valleys 
below the onset of stripes, but only displayed stripe pat- 
terns above T^^^ = 2.15 ±0.01, roughly 10% higher than 
in continuum simulations That study hypothesized 
that this discrepancy may be due to fluctuations which 
were unaccounted for in the continuum model. 

In Rayleigh-Benard convection of fluids near the onset 
of convection patterns, fluctuations caused by thermal 
noise create deviations from the dynamics predicted by 
Navier-Stokes equations without a noise source. Fluctu- 
ating hydrodynamics (FHD) theory models these fluctu- 
ations by adding noise terms to the Navier-Stokes equa- 
tions [20-[22|- FHD theory accurately describes the dy- 
namics of fluids near convection onset [2^ [l^l . Experi- 
ments indicate that fluctuations due to individual grain 
movement play a larger role in granular media than do 
thermal fluctuations in ordinary fluids [25j . 

Extending FHD theory to granular media is not trivial. 
The noise terms derived by Landau and Lifshitz f20j treat 
fluctuations near equilibrium which are small compared 
to the hydrodynamic fields and do not provide for local 
energy loss due to particle inelasticity. However, granular 
shaker experiments show fluctuations much larger than 
in ordinary fluids j25| , and any fluidized granular system 
is far from equilibrium due to inelastic particle collisions. 
In the shaken layers considered here, the mean free path 
of a particle is on the order of a particle diameter or less, 
so fluctuations due to small number statistics may be sig- 
nificant. Finally, recent simulations of a dilute granular 
gas [2^ showed that Landau-Lifshitz theory underesti- 
mates fluctuations in a ID homogeneous cooling state by 
neglecting memory effects of inelastic particles. 

As a test of the applicability of FHD, I treat fluctua- 
tions in the granular system analogously to thermal fluc- 
tuations in ordinary fluids. At each timestep, the sim- 
ulation calculates random local stresses and heat fluxes 
given by Landau and Lifshitz [20] at each grid point with 
no fit parameters, and includes these terms in the con- 
tinuum equations [13] . 

To visualize peaks and valleys formed by standing wave 
patterns, I calculate the height of the center of mass of the 
layer, Zcm (x, y, t) as a function of horizontal location in 
the cell at various times At a given time to and horizon- 
tal location {xi^.yi^')^ z^m 

{xo,yo, to) is the center of mass 
of all particles whose horizontal coordinates lie within a 
bin of size Axbin x ^yun centered at (a;o,yo)- The sim- 
ulation grid size defines the bins: /S.Xhin — Ayun — 2cr. 
Throughout this report, I characterize the patterns at the 
beginning of a cycle, when the plate is at its equilibrium 
position and moving upwards. Peaks in the pattern cor- 
respond to maxima of Zcm] valleys correspond to minima. 

An example standing wave stripe pattern is shown 
in Fig. [1] Continuum simulations both with (Fig. Ilbp 
and without noise (Fig. llap produce stripe patterns for 
r = 2.2 and /* = 0.4174. These patterns oscillate sub- 
harmonically, repeating every 2//, so the location of a 
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FIG. 1: An overhead view of a layer of grains, showing the 
center of mass height Zcm as a function of horizontal position 
[x, y) in a cell with horizontal dimensions L^xLy = 42a" x 42a, 
from continuum simulations (a) without noise and (b) with 
noise. Peaks of the layer corresponding to large center of mass 
height Zcm are shown in white; valleys corresponding to low 
Zcm are shown in black. 

peak of the pattern becomes a valley after one cycle of 
the plate, and vice versa [13]. When the accelerational 
amplitude is reduced to P = 1.9, stripes do not appear. 

In both cases, two wavelengths fit in the box for this 
box size and frequency (Fig. [TJ, although simulations 
without noise show sharper peaks and valleys with larger 
amplitude than simulations with noise. To compare these 
amplitudes, I examine the deviation of the height of the 
center of mass of the layer as a function of horizontal lo- 
cation in the cell from the center of mass height averaged 
over the entire cell: 

i^ix, y, t) = Zcm{x, y, t) - {zcm{x, y, t)) , (1) 

where brackets represent an average over all horizontal 
locations in the cell at a given time t. Thus, ^■0^(i)) 
represents the mean square deviation of the height of the 
layer from the mean height of the layer. Note that ("0^) 
is large for layers with high amplitude peaks and valleys, 
and goes to zero as the layer becomes perfectly flat. 

To distinguish between ordered patterns (stripes) and 
disordered fluctuations, I characterize the long range or- 
der of the pattern. I first calculate the power spectrum 
of the pattern S {kx, ky, t) as a function of wavenumbers 
kx and ky. Transforming to polar coordinates kr and kg 
in k space and integrating radially yields the angular ori- 
entation of the power spectrum S{kg). I bin kg into 21 
bins between kg = and kg — tt, and characterize the 
long range order by the fraction of the total integrated 
power that lies in the bin with the maximum power: 

p _ S„iax /n\ 

Jo ^ (kgjdkg 

where Smax is the integrated power within an angle 7r/21 
of the maximum value of S{kg). For a perfectly dis- 
ordered state, with equal power in all directions, Pmax 
would approach ^ ft! 0.05, while Pmax = 1 for a state 
with all power in a single bin. Thus Pmax provides a 
measure of order when stripes form. 
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FIG. 2: The mean square deviation {tp^) of the local center 
of mass height from the average center of mass height of the 
layer as a function of accelerational amplitude F for simu- 
lations without noise (a), and with noise (b). In (a), {ip^) 
is averaged over 50 cycles of a single simulation for each F 
(dots). The shaded region 1.94 < F < 1.96 indicates the 
transition between flat layers and layers with non-negligible 
peaks and valleys. For simulations with fluctuations, the data 
in (b) are averages (dots) with root mean square deviation 
(bars) from 50 cycles from each of six trials within the range 
1.90 < F < 2.20, and each of three trials outside that range. 



FIG. 3: Global ordering Pmax as a function of nondimensional 
accelerational amplitude F for continuum simulations without 
noise (a), and with noise (b). For simulations without noise, 
Pmax is averaged over 50 cycles from a single simulation and 
represented as dots, while for simulations with noise, Pmax 
is averaged (dots) over multiple simulations, with error bars 
calculated as root mean square deviation from this average. 
In both cases, there is a transition (shown in gray) to an ap- 
proximately constant Pmax ~ 0.4. The transition area shown 
in gray is 1.94 < F < 1.96 in (a) and 2.07 < F < 2.17 in (b). 



I examine (V'^) ^nd Pmax for simulations vifith varying 
F. In each case, the simulation begins with a flat layer 
above the plate with small amplitude initial random fluc- 
tuations. The simulation runs for 400 cycles of the plate 
to reach a periodic steady state. Then ("0^) ^-^d Pmax are 
averaged over the next 50 cycles. Compared to simula- 
tions without noise, simulations with noise show greater 
variation between cycles in their final state; I run these 
simulations three times for each F to find an average less 
influenced by transient behavior. As patterns occur for 
r = 2.20, but not for F = 1.90, three additional simula- 
tions (for a total of six) were run for each F in the range 
1.90 < F < 2.20 to more precisely locate pattern onset. 

For simulations without noise, fluctuations in the ini- 
tial condition decay over time for F < 1.94, producing a 
flat layer (Fig. [2a|) . As F increases, there is a jump to a 
periodic state of non-neglig ible (V-^) for F = 1.96, and 
large amplitude waves occur for all F > 1.96 (the region 
1.94 < F < 1.96 is shaded in Fig. [2al). When noise is 
added, the layer remains flat for some values F > 1.96 
(Fig. I2bl) . Non-negligible amplitudes of (tp^) occur for 
F > 2.0, but there is not a sharp jump in amplitude. 

bmce {^^) in Fig. [2b] increases gradually with increas- 
ing F rather than showing a sharp onset of waves, I ex- 
amine the order parameter Pmax to distinguish between 
stripes and disordered fluctuations as shown in Fig. [31 
For simulations without noise, all layers with F > 1.96 
show a nearly constant value of Pmax ~ 0.4 (Fig. [3a| . 
corresponding to the stripe patterns seen in Fig. [Taj For 
F < 1.94, the initial fluctuations decrease over time, lead- 
ing to a very flat layer (c/ Fig. [2a|) with lower Pmax- I 
identify the critical value Fg*"* = 1.95±0.01 above which 
stripe patterns are formed in simulations without noise. 

For noisy simulations, there is relatively large uncer- 
tainty in Pmax in the shaded region 2.07 ^ T < 2.17 
(Fig. I3bl) . Visual inspection shows transient behavior in 



this region, with temporary order appearing and then 
disappearing, yielding variation in Pmax from simulation 
to simulation. Above this shaded region, Pmax ~ 0.4 
with low variation, indicating consistently reproducible 
stripes. Below this region, Pmax is consistently lower, in- 
dicating disordered fluctuations. I thus identify the criti- 
cal value above which stripe patterns form in simulations 
with FHD terms T^"^ = 2.12 ± 0.05. 

These results for continuum simulations without noise 
pcont _ 2 95 jj- Q Qi agree with previous continuum sim- 
ulations showing an abrupt transition from a flat layer 
to stripe patterns at Fg""* = 1.955 ± 0.005 Simula- 
tions with FHD noise, however, show a gradual increase 
of disordered fluctuations below the onset of ordered 
stripes, and a transition to stripes at T^^^ — 2.12±0.05. 
While continuum simulations with noise differ from those 
without noise, they are consistent with previous MD 
simulations showing the transition to stripe patterns at 
f"cf ^ — 2.15 ± 0.01, with a gradual increase in amplitude 
of disordered fluctuations below this value [4]. 

Finally, I investigate the wavelengths of these patterns. 
Experiments have shown that wavelength A depends on 
the frequency of oscillation [2^ [2^ . For a range of layer 
depths and oscillation frequencies, experimental data for 
frictional particles near pattern onset were fit by the func- 
tion A* = 1.0-h l.l/*-i-^2±o.03^ y ^ x/H [H. 

I investigate frequency dependence by holding dimen- 
sionless accelerational amplitude F = 2.2 constant, while 
varying dimensionless frequency /*. Simulations were 
conducted in a box of size = 168(t, Ly — 10a, and 
Lz = 160(7. This orientation causes stripes to form par- 
allel to the y— axis. The dominant wavelength was calcu- 
lated from the wavenumber kx in the x— direction which 
exhibited the maximum power during 50 cycles of the 
oscillatory state. Due to the periodic boundary condi- 
tions and finite box size, wavelengths must fit in the box 
an integer number of times, yielding uncertainty in the 
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FIG. 4: Dispersion relations for stripes which form perpendic- 
ular to the long dimension of cells with horizontal dimensions 
168(T X lOcr. Data for simulations with noise are shown as 
squares, without noise as triangles, and points where the two 
simulations yield the same wavelength are shown as circles. 
Error bars are calculated exclusively from discretization due 
to periodic boundary conditions in a finite size box. In both 
simulations, the dominant wavelength of the final oscillatory 
state A fits very well to the dispersion relation found in ex- 
periments A* = f.0 + f.l/"^-^^±°"^ (solid line) 

wavelength that would be selected in an infinite box. 

For this box size, frictionless MD simulations and con- 
tinuum simulations without noise have been shown to 
agree with experimental results for frictional particles 
through the range 0.20 ^ ft < 0.45; friction appears 



unimportant in wavelength selection through this range 
Wavelengths found in continuum simulations with 
and without noise are compared to the dispersion rela- 
tion fit to experimental data in Fig. 2] Both simulations 
agree quite well with the experimental fit throughout this 
range. The addition of noisy fluctuations does not appear 
to significantly affect the wavelength of the patterns. 

In conclusion, continuum simulations without friction 
can describe important aspects of pattern formation in 
granular media. With or without noise, frictionless con- 
tinuum simulations produce patterns with wavelengths 
consistent with experimental results in layers of particles 
with friction. For the shaken layers studied in this report, 
patterns in continuum simulations without noise occur 
for critical accelerational amplitude Tc approximately 
10% lower than in experimentally verified molecular dy- 
namics simulations. Including fiuctuating hydrodynam- 
ics (FHD) alters the onset of patterns; Tc for continuum 
simulations with noise is consistent with MD simulations, 
but not with continuum simulations lacking this noise. 
These results indicate that fiuctuations play a significant 
role in this system, and also suggest directions for fur- 
ther research. Simulations including memory effects 
or other variations in the FHD model could be compared 
to test which approximations significantly alter pattern 
formation. In addition, testing the effects of noise on 
other granular systems will be important in establishing 
a general theory of granular hydrodynamics and the role 
of fluctuations within that theory. 
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